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Motivated by observations of the dynamics of Myxococcus xanthus, we present a self-interacting 
random walk model that describes the competition between chemokinesis and chemotaxis. Cells 
are constrained to move in one dimension, but release a chemical chemoattractant at a steady 
state. The bacteria senses the chemical that it produces. The probability of direction reversals is 
modeled as a function of both the absolute level of chemoattractant sensed directly under each cell 
as well as the gradient sensed across the length of the cell. If the chemical does not degrade or 
diffuse rapidly, the one dimensional trajectory depends on the entire past history of the trajectory. 
We derive the corresponding Fokker-Planck equations, use an iterative mean field approach that we 
solve numerically for short times, and perform extensive Monte-Carlo simulations of the model. Cell 
positional distributions and the associated moments are computed in this feedback system. Average 
drift and mean squared displacements are found. Crossover behavior among different diffusion 
regimes are found. 

PACS numbers: 87.17-d, 87.17- Jj 



CO 



I 

C 

o 
o 



> 

00 

m 

OS 

o 
m 
o 



-a 
c 

o 
o 



X 



I. INTRODUCTION 



The dynamics and pattern formation of bacteria serve 
as paradigms in understanding many properties of multi- 
cellular interacting systems, such as collective behavior, 
self organization, evolution, and development 0,0- The 
characteristics of bacterial motility and aggregation de- 
pend on numerous biological and chemical parameters, 
which include light exposure, temperature, concentration 
of food or of other substances, and lead to the existence 
of many classes of cell motion. Since a bacterium can- 
not utilize electromagnetic or acoustic radiation to sense 
its environment, it must rely on physical contact and its 
motility may depend on the production, diffusion and 
degradation of a relevant set of sensed chemicals. These 
chemicals are called chemoattractants or chemorepellants 
if they attract or repel bacteria, respectively. Exam- 
ples of chemoattractants include food - sugars and amino 
acids - whereas antiobiotics, fatty acids or other noxious 
substances are chemorepellant. 

A prototype and well studied example of bacterial 
motility is the run and tumble mechanism used by Es- 
cherichia coli in liquid environments. The specialized 
structures allowing for cellular motion are known as flag- 
ellae, spinning helical tails which extend from the cellular 
membrane into the surrounding environment: counter- 
clockwise rotation of the flagellar motor leads to a one- 
directional run whereas clockwise rotation causes tum- 
bling along a random direction. During its run, an E. coli 
cell periodically senses for chemosensitive substances, 
and, by comparing concentration levels in new and old 
environments, adjusts its motion and its likelihood to 
tumble in a new, randomly chosen, direction 0. 0]. If 
the cell is moving in the direction of increasing attrac- 
tant, for instance, the probability of tumbling is lowered, 
so that the run time and 'mean free path' length in this 
direction are increased. An E.coli cell, thus, compares 



chemical concentrations at nearby points, so that its fu- 
ture motion is determined by a chemical gradient 0, 
via the mechanism of chemotaxis. 

Chemotaxis may occur in other non-acqueous systems, 
such as in colonies of the bacterium Myxococcus xan- 
thus or of the eukaryote Dictyostelium discoideum (slime- 
mould) which crawl on agar plates @|. The exact regu- 
latory mechanism leading to the surface gliding of these 
rod-like and slow moving organisms is yet unknown, and 
responses to external stimuli appear to be more com- 
plex than in enteric bacteria such as E.coli 8J. The 
cells undergo chemotaxis, but lack, or show very mod- 
est, responses to specific nutrient stimuli 0, Hcf . sug- 
gesting that the chemotactic behavior is due mainly to 
self-generated signaling chemicals [ll| . The fact that cell 
density regulates typical reversal frequencies may be a 
consequence of one bacteria sensing the higher levels of 
signaling chemical in the presence of others |12| . In par- 
ticular, under starving conditions, the cells are driven by 
auto-released fibril trails to aggregate and form multi- 
cellular compounds called fruiting bodies. These com- 
pounds raise above the agar plate and eventually sporu- 
late into new and more resistant organisms which are 
then released to the environment in search of new sources 
of nutrients. 

One major difference between these systems and en- 
teric bacteria is that the chemoattractants embedded in 
or self-deposited on the agar plates are relatively im- 
mobile on the time scale of bacterial motility, whereas 
chemoattractant substances in acqueous environments 
have a finite diffusion constant. Like E.coli, the gliding 
cells can only sense what they come into immediate con- 
tact with, and since the chemosensitive substance diffuses 
slowly, large fibril concentration differences exist over the 
length of the cell. The gradient driving the chemotaxis 
mechanism is then determined between positions within 
a cell body length, as opposed to the large distances in- 
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volved in enteric bacteria chemotaxis. 

Another possible mechanism driving cellular motility 
is chemokinesis. Here, changes in the direction of mo- 
tion or in the cell velocity are determined locally, with- 
out recourse to the determination of concentration gra- 
dients. Chemokinesis may be a function of temperature, 
substrate adhesion, salt concentration or substances af- 
fecting the internal metabolism of the cell, slowing or 
enhancing its motion or turning frequency . 

The observed patterns of movement of M. xanthus and 
of slime-mould suggest that the motility of these systems 
may be affected by self-induced chemokinesis as well as 
self-induced chemotaxis. M. xanthus move by shooting 
out and retracting pili from the ends of their prolate bod- 
ies |3 113 • I n nutritionally abundant environments or 
under conditions of high population density, both cellu- 
lar systems generate repulsive fibril trails, leading to cell 
dispersion and the colonization of new regions of the agar 
plate. The mechanism is hypothesized to be dictated by 
chemokinesis [l6[ll7j. 

Previous theoretical studies into the collective dynam- 
ics of interacting motile bacteria, both of E. coli and M. 
xanthus type, include coarse-grained convective-diffusive 
transport type models flM ITa. I20L I2TL |22| and lattice- 
type simulations [2^, 124 125| . These studies have consid- 
ered the effects of interactions with concentration fields, 
modeling chemotaxis or chemokinesis or, as in Ref. |26| . 
both chemosensitive mechanisms. The continuum theo- 
ries take into account turning rates, and changes in speed 
or direction of motion, however they generally ignore long 
time path history effects induced by the sensed chemicals. 
In these systems, the bias to the motion arises from an 
external, fixed stimulus and not a dynamical one. Similar 
studies - with a fixed external trail - have been presented 
in different contexts, such as the aggregation of ants due 
to chemotaxis J27|. Other authors [23 have presented 
systems in which the fibril trail that dictates the motion 
of the cells is determined over an extended area and is 
not sensed locally. 

The aim of this paper is to model and disentangle 
the effects of self- generated chemotaxis and chemokine- 
sis, acting individually or concurrently and under diverse 
attractrant and repulsive conditions. In contrast to many 
of the previous studies, the chemical trail of our systems 
will emerge from the motion of the bacterial cells and 
evolve with the cellular motion itself. Observations of 
individual organism paths may be more revealing of cel- 
lular dynamics than multicellular ensembles, and we re- 
strict our attention to the behavior of a single bacterium. 
Even if isolated, the dynamics of this single cell will be 
affected by its self-secreted fibril trail, and the resulting 
pattern will be that of a self-interacting random walker. 
We shall try to understand cellular motion by developing, 
and analytically solving, the equations of motion for an 
isolated M. xanthus cell under a mean-field approxima- 
tion, and by means of simple Monte-Carlo simulations, 
incorporating chemotaxis and chemokinesis to various de- 
grees. The underlying assumption of this paper is that 
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FIG. 1: The trajectory of a single bacterium of width w and 
length i crawling in one dimension and depositing chemoat- 
tractant (or repellent) on the surface. The density of shading 
of the footprint is proportional to the total attractant de- 
posited and arises from the hypothetical trajectory depicted 
by the arrows beneath. For the sake of simplicity, we have 
neglected the cell length £ in the depiction of the trajectory. 



the chemosensitive material acts as a regulator of cell mo- 
tion affecting only the frequency of direction reversal, via 
chemokinesis or chemotaxis. The speed of the bacterium 
is assumed constant in either direction and is controlled 
by inherent cellular bio-energetics. 

Experimental single cell trajectories show that an iso- 
lated cell travels along its axis, occasionally veering off, 
and for the sake of simplicity, we will only consider one- 
dimensional dynamics. In order to validate our models, 
cell particle tracking can be eventually performed in one- 
dimensional etchings or imprints into the agar plates. 



II. CONSTANT REVERSAL RATES 

In this Section we derive a general Fokker-Planck equa- 
tion for the single one-dimensional random walker. As 
discussed in the Introduction, we assume that the 'bac- 
terium', or particle, travels with constant speed i>o, and 
is subject to an ad hoc directional reverse mechanism 
that incorporates both chemotaxis and chemokinesis. To 
start, we ignore the effect of both phenomena, and we 
simply investigate the role of a constant reversal rate. 

The probability that a bacterium is centered at posi- 
tion x at time t with velocity ±w will be denoted by 
P±(x,t) and the probability of a directional switch in 
time dt is dt/r — 7 dt. If 7 is independent of (x,t) we 
can write: 

P+(x,t) = P+(x-v dt,t-dt){l-~f + dt)+ (1) 

P-(x,t)^dt, 
P-(x,t) = P_(x + v dt,t-dt)(l-j_dt)+ (2) 

P+(x,t) 1+ dt. 

In the above equations, we have included the possibility 
that 7 depends on the type of reversal, since in general 
the reversal rates from one direction to the other are not 
the same. Specifically: 

7+ : +v -> -v , 7_ : -v -> +v - (3) 
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Expanding Eqns. and Q and keeping only the 0(dt) 
terms, we find: 

P+(x,t)+v d x P+(x,t) = - 7+ P+(x,t)+7_P_(x,t),(4) 
P-(x,t)-v d x P-(x,t) = - 1 -P-(x,t)+j+P+(x,t).(5) 

The initial conditions are chosen so that at t = the par- 
ticle is localized at x = with amplitudes a± satisfying 
the condition that a + + a_ = 1: 

P+(x,0) = a+8(x), P-(x,0) =a-8(x), (6) 

We obtain: 

(q-x - q+v t) 



P+(x,t) 



exp 



2w 



(7) 
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Here, the parameters g+, q- and Q are defined so that 
q+ = (7++7_), q- = (7- -7+) and (q+Q) 2 = (^-g 2 ). 
The function O is the Heaviside function, 9(.t) = 1 if 
x > 0, and the / functions are modified Bessel func- 
tions. The 6 function terms carry the probability of the 
random walker cell having performed no reversals, and 
the Heaviside terms embody causality. 

The total particle probability distribution function is 
P(x,t) = P + (x,t) +P_(x,i). Under the symmetric con- 
dition a+ = a_ = 1/2, and 7+ = 7_ (or equivalently 
<7_ = 0) the solution reduces to that for the generalized 
Smoluchowski diffusion equation determined by Hemmcr 
|29|. which at large times is a spreading Gaussian. The 
Fourier transform of the distribution functions yields the 
moments of the random walker cell from which the cu- 
mulants can be extracted. 



III. CHEMOKINESIS AND CHEMOTAXIS 

We now consider the generation and sensing of 
chemoattactant or chemorepellant, and its effect on the 
reversal rates 7±. M. xanthus cells deposit attractant 



chemical matter contained within a fibril slime under the 
area of cell contact. We assume that the release of the 
chemoattractant occurs at a constant rate and uniformly 
over the cell-substrate footprint. In general, for a cell of 
length £, the fibril concentration 4>(x, t) obeys the follow- 
ing equation: 

<P(x,t) = DV 2 c/)(x,t)+ (9) 
rG[I/2-|a-X(t)|]-r^(x,t), 

where the chemical is generated at a constant rate r per- 
unit area. The Heaviside function in Eq. © limits 
the increase in fibril concentration to the region strictly 
under the footprint of the bacterium centered at X(t), 
whereas degradation and diffusion are included via the 
terms proportional to and D. If the surface-deposited 
chemoattractant does not diffuse appreciably and PV 2 
is neglected, we find: 



(x,t) = <j>{x,Q)e~ rdt + (10) 
r[ dt'e r « (t '- t '>e[l/2-\x-X(t')\}. 



Note that <f>{x, t) depends on the entire history of the 
trajectory X(t') up to time t' = t. To couple the fibril 
concentration with cell dynamics, we must now include 
a relation between "f±(x, t) and <f>(x, t). This relationship 
will depend on the particular mechanism of cell direction 
reversal and various choices will be discussed in the fol- 
lowing subsections. It is important to note that, since 
we are assuming that fibril sensing is the main factor 
in determining direction reversals, we are also implicitly 
assuming that the sensing process, for example via a bio- 
chemical network within the bacterium, is much faster 
than the typical time for a cell to have moved apprecia- 
bly. If we further take the molecular sensing apparatus to 
be distributed uniformly at the cell substrate footprint, 
the total chemoattractant <&(X(t\t) sensed at time t is 
just the local contribution integrated over the cellular 
length: 



X(t)+£/2 



c/)(x,t)dx, (11) 

ix(t)-e/2 

where w is the width of the cell as depicted in Fig.Q] 

A. Chemokinesis 

Let us consider the case of chemokinesis first. Here, 
the reversal rate depends only on the integrated fibril 
concentration. We will assume no intrinsic drift and 
set j{x,t) = ■y + (x,t) = 7_(ar, t). It is also reasonable 
to expect the switching probability to saturate for large 
enough fibril concentration, when the cell cannot respond 
to further increases in fibril levels. A plausible functional 
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dependence for 7 on the total sensed chemoattractant $ 
is of the type found in cooperative chemical binding, a 
Michaelis-Menten form with Hill coefficient /3, so that for 
a cell at position X(t) = x, the reversal rate is: 



7(35, t) = 70 + K(x,t), 

= 70 + tip- 



in) 



where $0 an d /?, the inflection parameter and transition 
sharpness (or Hill coefficient), are peculiar to the system 
being considered. The effect of chemokinesis is measured 
by 8p; with the unbiased choice 70 = 1/2, Sp is limited 
to -1/2 < 5p < 1/2. 

We must now solve Eqns. (0J and (J5J with the above 
spatial and temporal dependence for 7(2;, t), which in 
turn depends on the concentration $(x,t) via Eqns. i|lU|) 
and 1|11[1 . The latter equations carry an explicit depen- 
dence on the past history of the walker. Our approach 
will be to consider the 70 contribution of Eq. l(T2*|> as giv- 
ing rise to Eqns. 10} and JSJ with constant reversal rates, 
and find its associated Green's function. The k(x, t) term 
in j(x, t) of Eq. (|12|) then generates non-homogeneous dif- 
ferential equations, whose solution may be obtained via 
the known Green's function. We will then be able to con- 
struct an iterative process to obtain P(x,t) as a function 
of all preceeding probabilities P(x' , t' < t). To build the 
Green's matrix function G(x,t), we note that Eqns. Q 
and JSJl can be written as: 



P(x,i) = M(a;) POM), 
where P(x,t) = [P + (x,t), P-(x,t)f and 



(13) 



M(x) 



-v d x - 7+ 7_ 
7+ v d x - 7_ 



(14) 



For the case of generic 7±, independent of (x,t), the so- 
lution to this equation, with the initial conditions ex- 
pressed by Eqns.® is given by Eqns. and (JHJ. In 
the case 7+ = 7_ = 70 we refer to these solutions as 
P{ 0) (x,t) and P W (x,t). The Green's function for this 
problem stems from the modified matrix equation: 



leads to the conclusion that the Green's matrix function 
is given by: 



G(x,t) 



Pi 0) (x,t) a+=1 Pl 0) (x,t) a+=0 \ 



This result can also be verified by direct evaluation of 
G(k,u>) in Fourier space. From Eqns. 0} and 10 it 
follows that reversal rates of the type r y±{x,t) = 70 + 
K±(x,t) yield a new matrix equation: 



P(x, t) = M(x)P{x, t) + K(x, t)P(x, t). (19) 



where: 



K(x,t) 



— K + {x,t) K-(x,t) 
K+{x,t) —K_(x,t) 



(20) 



The resulting equation for P(x,t) can be solved itera- 
tively: 

P{x,t) = pW(x,t)+ (21) 
I dt' dx' G(x-x',t-t') K{x',t') P(x',t'). 

JO J ~oc 

In the case of chemokinesis, the non-homogeneous term 
K+(x,t) = K,-{x,t) is contained in Ea. I|12|l . In order 
to simplify the expression for <i>(x,t), we assume unit 
cell width and take the limit of cellular length I — > 0. 
We also neglect diffusion and decay of the fibril. Under 
these assumptions, 4>{x,t) = r£S(x — X(t)) and Q(x,t) = 
w£(j)(x,t). The problem is now a deterministic one, as 
exemplified by the delta-function in the expression for 
cf>. In order to apply Eq. (|21|l . which gives the statistical 
probability for one bacterium to lie at (x,t), we need to 
know its exact location at X(t') for all previous times. In 
other words, the complete history of the cell is needed to 
determine future motion and solutions cannot be found. 

Nevertheless, we can approximate P±(x, t) by an aver- 
aged density p± , utilizing a mean field theory approach. 
Here, we average the probabilities P± for a single walker 
over many distinct independent realizations, so that the 
self interaction is to be expressed on the average of all 
replicas of the system, and does not depend on the his- 
tory of individual walkers. We can now rewrite the equa- 
tion for the fibril concentration as: 



G(x-ar',t>0) = M(x) G(x - x',t) + 

18(x-x')5(t), (15) 
G(x-x',t< 0) = 0. (16) 

The fact that: 

P (0) M) = ( ^o) J M = G(x, t) ( °+ ) , (17) 



<fi(x,t)=rtp(x,t), (22) 

where we indicate the total bacterial density by p(x,t). 
The above relationship signifies that, on average, fibril 
growth is proportional to the total density of bacteria at 
(x,t). Integrating Eq. 122|) and taking the initial condi- 
tion 4>(x, 0) = 0, we obtain: 

<j>{x,t) =r£ f dt'p{x,t'), (23) 
Jo 
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from which, using Eq. I|21(l and the definition Ap(x, t) 
p + (x,t) — p_(x,t) we can write: 



systems. In the limit of cellular length t 
expressions are written as: 



0, the above 



Ap{x,t) = Ap (0 \x,t) + 



Sp / dt' 
Jo 



dx g(x — x ,t — t) 



(24) 

r Ap(x',t')^(x',t') 



$ + 4>P(x',t>) 



Here g(x, t) is a combination of the matrix elements of the 
Green's function: g ~ —gu + 9n + .921 — 322 and is even 
in x. An equation similar to Eq. (|24(l can be obtained for 
the total density p(x,t), which depends on Ap(x,t): 



p(x,t) = p (0 \x 1 t) + 

5p / dt' / dx'h(x -x',t- t') 



(25) 

Ap(x',t')^(x',t') 
4 + ^{x',t') 



Here, h — —gu + gu — 321 + 322 and is odd in x, as re- 
quired by the normalization of both p{x,t) and p°(x,t). 
It is useful to note that the above recursive equations 
may be used to evaluate the number density n(x, t) for a 
system of N interacting bacteria as well. In this case, the 
distribution function n±(x,t) = Np±(x,t) and the evo- 
lution equations are the same as l|24|l and (|25|l provided 
we redefine O = Nn . Equations (|23|) . 1241) and (|25|l can 
now be solved numerically. Once Ap{x, t) and p(x, t) are 
determined, <fi(x, t) can be evaluated and used to calcu- 
late the quantities of interest at subsequent times. 



B. Chemotaxis 

The incorporation of chemotaxis into the determina- 
tion of the reversal rate follows closely that of chemoki- 
nesis. Let us consider chemotaxis as the sole reversal 
mechanism for a cell of length I traveling at vq speed in 
the positive direction, and situated at position X{t) = x. 
The chemical gradient over the cellular length will drive 
the reversal rate. To specify the relation between gra- 
dient and reversal probabilities, we first introduce the 
probability terms p±(x,t): 



p±{x,t) = exp [=fcr<f> x (x,t)] 



(27) 



where the 4> x = dcj)/dx. We use this exponential function 
to capture the sensitive chemical signalling, such as that 
found in E. coli 30] . Another alternative for the function 
p± is to include a saturation in the form of <j) x /4> in the 
exponent of (|27|l . 

We assume the reversal probabilities r y±(x,t) at posi- 
tion (x,t) to depend on p±{x 1 t) through the following: 



7±0M) = 2 



fi±(x,t) - 1 



p±(x,t) 



(28) 



the 



so that for a = 0, in the absence of chemotaxis, 
reversal probabilities are 7± = 7o = ^ • 

We are now able to use the same formulation de- 
rived for the chemokinesis mechanism and apply it to the 
chemotactically driven case. Again, we must resort to the 
mean field calculation of the densities p(x, t) and Ap(x, t) 
by means of the Green's matrix function. The matrix K 
appearing in the non-homogeneous term of Eq. (|19|l now 
has components K±(x,t) given by: 



K±(x,t) 



p±(x,t) - 1 
p±{x,t) + 1' 



(29) 



and Eq. I|21|l still holds. Note that since p-(x,t) — 
p^i^Xjt), then n + (x,t) + K_(x,t) = 0. Repeating the 
same type of calculations as for the chemokinesis case, 
for Ap under chemotaxis we find: 



Ap(x,t) = V°W) + 2 



(30) 



dt' / dx'g(x-x',t-t')p(x',t')K + (x',t'). 



p±(x, t) = exp 



(j){x +£/2,t) - <f>(x - e/2,t) 



(26) 



where a is a constant that measures the strength of 
chemotaxis. For positive cr, the above expression trans- 
lates into a low reversal probability whenever cj){x + 
1/2, t) > (f>{x - 1/2, t) so that the cell is likely to keep 
moving speed at +vo. This means we are modeling an 
attractant fibril trail. On the other hand, negative values 
of a will represent a chemorepellant system. The same 
argument can be presented for a bacterium traveling at 
— v speed. In this case, the direction of motion tends 
to persist for negative gradients, hence the expression for 
p~(x, t), where positive a values signify chemoattractant 



Here the previous history of the cell is contained in 
k + (x' , t') through the derivatives of the fibril concentra- 
tion cj)(x',t'). The g function is the same as defined in 
the case of chemokinesis. Similarly, the density p(x, t) is: 



p(x,t) = p {Q) (x,t) + 



1 



(31) 



dt' / dx'h(x -x',t- t')p{x', t')K+(x', t') 
J-00 



and p(x, t) is properly normalized by the fact that h(x, t) 
(defined in the preivous subsection) is odd with respect 
to x. 
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C. Combined Chemokinesis and Chemotaxis 

If the two mechanisms of chemokinesis and chemotaxis 
are active, both bp and a are non-zero, and j±(x, t) must 
include both contributions. For 7±(x,t), we write: 



7±0M) = g i 1 



Pk(x,t)[n±(x,t) + 1] - 1 
Pk(x,t)\p,±(x,t) - 1] + 1 



(32) 



Here, pk(x,t) is the contribution from chemokinesis in 
the form expressed by Eg. H12[) with 70 = 1/2, i.e. 
Pk{x, t) = 7(2;, t) from the chemokinesis case. The contri- 
bution from chemotaxis is fi±(x, t) as defined by Eq. (|2"?jl . 

The choice a — sets 7±(x,i) = pk{x,t), that is, 
the rate is dictated only by the chemokinesis mechanism. 
Conversely, 5p = sets pk = 1/2 and j±(x,t) reduces to 
the purely chemotactic expressions of Eq. 128|) . For both 
(j = Sp = the reversal rate is 1/2. The corresponding 
n + (x,t) and K-(x,t) are contained in Eq. (|32|l and yield 
the following recursive relationship for Ap(x,t): 



Ap(x,t) = Ap {0) (x,t) + / dt' dx'g(x-x',t-t') 

JO J — oo 

\ [U(x', t')p{x\t') + V(x', t')Ap{x>, t')] . (33) 
The corresponding equation for p(x, t) is: 



p(x,t) = p {Q) (x,t) + / dt' / dx'h(x-x',t-t') 

JO J -oc 

~ [U(x', t')p(x', + V(x', t')Ap(x', t')} . (34) 

Here, the U{x, t) and V(x, t) functions depend on chemo- 
taxis and chemokinesis via: 



U 



V 



[p k (^+ - 1) + (Pfe - 1) -Pfc]' 

~ 2p fc ) 

[Pfc(M+ - 1) + - !) -Pfc] ' 



(35) 



(36) 



We have used /i+^t- = 1 and suppressed the {x,t) depen- 
dence. Note that in the absence of chemotaxis, for a = 
and p + = 1, U = and V = 2p fc — 1, leading to the pure 
chemokinesis case. In the absence of chemokinesis, for 
Pk = 1/2, IA = k + , as defined in the previous subsection, 
and V = 0. We shall discuss the numerical solution to 
the equations for p(x, i) and Ap{x, t) in the various cases 
of chemokinesis, chemotaxis or both, in the next Section. 



IV. DISTRIBUTION FUNCTIONS 

In this Section, by means of numerical integration, we 
solve equations (J23J, ED, (EOJ, ESI, and (O for 



chemokinesis 

4>„=0.2 p=3 



0.06 



o 

CO 



0.04 



0.02 



0.00 




FIG. 2: Distribution functions in the case of chemokinesis for 
4>o — 0.02 and P — 3. Note that Sp > increases the reversal 
rate, localizing the cell at the initial position and Sp < 
values have the opposite effect. The time scale is n = 600 
time steps in units of At = 0.05. 
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FIG. 3: Distribution functions in the case of chemotaxis 
for different a values. Note that a > values represent 
chemoattraction and the cells are localized at the initial posi- 
tion, a < values instead represent chemorepulsion and two 
reprelling opposite peaks develop. The time scale is n — 600 
time steps in units of At = 0.05. 



p(x, t) and Ap(x, t). Under chemokinesis, Eqns. (|25|l and 
(|24|l are Volterra equations of the second kind in the t 
variable, for which standard methods can be applied |3l| . 
In particular, we discretize both the temporal and spatial 
axis according to a uniform mesh, of spacing At and Ax, 
and solve the two coupled equations iteratively. 

Let us assume that Ap(x',t' < t) and p(x',t' < t) are 
known, where t' represent time steps up to t — At. At 
the subsequent time step t, Ap(x, t) is evaluated by cal- 
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FIG. 4: Fibril concentration 4>{x,t) in the case of either 
chemotaxis or chemokinesis for different a and Sp values. The 
time scale is n — 600 time steps in units of At = 0.05. In the 
chemokinesis curves cf>o = 0.2, close to the saturation limit of 
the Michaelis-Menten form of Eg .1121 



culating Ap^(x, t) and by adding to it the sum over all 
previous time steps t' < t of the integrand. Of these 
terms, the last one, at t — t', contains the kernel func- 
tion g(x — x',0) — —2S(x — x'), which, when integrated 
over space, extracts a term proportional to the quan- 
tity of interest, Ap(x,t). The other contributions of the 
sum, for t ^ t' ', are spatial integrals, that are evalu- 
ated as sums after the discretization over the x' variable. 
These sums contain the known values of Ap(x',t' < t), 
as well as k(x' , t') which depends on the fibril concentra- 
tion <fi(x',t'). The latter is just the integrated p(x',t') 
as follows from Eq. i|23|) . In order to evaluate 4>(x' , t'), a 
time discrctizcd summation over the p(x' . t' < t) values is 
carried out. The assumption here is that the bacterium 
senses the environment before laying down new fibril, so 
that the concentration <fi(x',t') is evaluated not by the 
trapeziodal rule, in which the extrema of the integration 
region are weighed each as 1/2, but by weighing t = 
with 1 and ignoring t = t' . In the case of chemokinesis, 
once Ap(x,t) is determined, p(x, t), which depends on 
Ap(x',t'), can be calculated in the same way. 

In Fig.[21we plot the results for p(x, t) for a system un- 
der chemokinesis. The meshes are set at Ax = At = 0.05 
and the maximum time step, in these units is t = 600. 
Cell velocity is fixed at Un = 1 and the parameter 
4>o = 0.2. For chemoattractant fibril, or 6p > 0, the 
distribution is narrower than the bare Sp = case, also 
plotted in the graph. Negative values of 5p represent 
chemorepulsion and broaden the distribution, eventually 
leading to a plateau-like region. The existence of the 
plateau is due to the fact that for large negative values 
of Sp ~ —1/2 the reversal rate at saturation is almost 
zero, leading to ballistic motion of the cell in the region 
where the fibril has been deposited. Once the cell reaches 



FIG. 5: Distribution functions under chemotaxis, set at a = 
+1.0, and chemokinesis set at various Sp values. The time 
scale is n — 600 time steps in units of At — 0.05. 



the edge of the plateau region of fibril concentration, the 
reversal rate increases dramatically and the cell tends to 
reverse its motion and travel backwards, until the oppo- 
site side of the plateau region is reached. This behavior 
accounts for the uniformity of the distribution function in 
the central region. Time evolution also tends to broaden 
- and consequently lower - the plateau region, since the 
bacteria tend to build up fibril at the edges, thus increas- 
ing the range of motion. 

For the case of chemotaxis, p(x, t), is easier to evaluate, 
since it can be calculated without recourse to Ap(x,t), 
as indicated by Eq. I|31l) . The procedure is the same as 
that outlined for the case of chemokinesis. We take the 
gradient appearing in the expressions for p± to be eval- 
uated over a distance I = 2, so that for x — x% and t = t - t 
of the mesh we have: 



p±(xi,tj) = cxp 



^[^{xi+Wtj) ~ ■ (37) 



Results are plotted in FigO As discussed in the previous 
section, positive (negative) values of a indicate chemoat- 
traction (chemorepulsion). It is interesting to note that 
for large negative values of a two peaks develop (only 
one is seen from Fig-El the other is its symmetrical im- 
age in the x axis), with the bacteria traveling towards 
unexplored regions of space, void of fibril. 

When both mechanisms are present, chemotaxis and 
chemokinesis may enforce or compete against each other 
in their localizing or spreading effects. Chemoattractant 
values of a > for instance, contrast the ballistic ten- 
dency of Sp < 0, and viceversa, chemorepulsive values 
of it < oppose the localizing effect of Sp > 0. In Fig- 
ures through {7\ we examine the interplay of the two 
mechanisms for several values of a and Sp. 
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FIG. 6: Distribution functions under chemotaxis, set at a — 
—0.5, and chemokinesis set at various Sp values. The time 
scale is n = 600 time steps in units of At — 0.05. 
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FIG. 7: Distribution functions under chemotaxis, set at a = 
—4.0, and chemokinesis set at various Sp values. The time 
scale is n — 600 time steps in units of At — 0.05. 



For both cr, Sp > the distribution functions are lo- 
calized about the central position, whereas for cr, Sp < 
two peaks traveling in opposite directions develop. The 
sharpness of these peaks increases with increasing |cr|, 
as can be seen by comparing Figures H3 through [7| and 
the localization position increases with increasing \Sp\ as 
each figure also shows. Competing effects arise when dp 
and a are of opposite signs. For large negative values of 
cr, the effect of positive Sp values is simply to shift the 
mean of the spreading peak towards the origin (Figures 
0) ■ For smaller negative values of a increasing Sp tends to 
localize the distribution functions at the origin, until the 



traveling peak is washed out and a plateau region, typical 
of large positive Sp values forms (Fig.[f)|). Finally, in the 
regime of cr > and Sp < 0, the localizing tendency due 
to chemoattraction is offset by the chemokinesis-platcau 
tendency. Unusual distribution function patterns may 
arise as in the case of cr = 0.5 or cr = 1 and Sp — —0.4 or 
Sp = -0.48 in Fig.Ha 



V. MONTE-CARLO SIMULATIONS 

To analyze the statistics of the one dimensional dy- 
namics of cell motion under chemotaxis and chemokine- 
sis, we also implemented a stochastic simulation on a ID 
lattice. We discretize space and time and assume the ran- 
dom walker cell to have finite length I. At time t = 0, the 
single cell is positioned at the origin of a one dimensional 
track and there is no fibril substance present. Fibril is 
constantly secreted by the cell and the dynamics of the 
particle position obeys: 



X t+ \ = X t + T) t (X t - X t -i), 



(38) 



where X t +\ is the position of the center of the cell at 
time t + 1 and X t — X t -\ indicates the direction it was 
traveling prior to time t. The stochastic random variable 
r]t is defined as: 



-1 with prob. jt 
-1 with prob. 1 — "ft 



(39) 



where "ft is the probability for the cell to reverse its di- 
rection, before the time step at t + 1. The initial direc- 
tion of motion is chosen so that X_i is ±1 with equal 
probability. We simultaneously model chemotaxis and 
chemokinesis through the dependence of jt on the self 
secreted fibril concentration (f>(x,t), as in the continuum 
case. Discretizing equations © and we obtain: 



<j> t {x) = 5^0(4/2-|a;-jr,|), and (40) 

8=0 

e/2 

$ t (x) = ]T M*t + x), (41) 

x=-£/2 

where 4>t{x) is the amount of chemoattractant that has 
built up at location x and $t is the total amount of 
chemoattractant sensed by the cell. Both 4>t(x) and $t 
naturally depend on the past trajectory of the cell, which 
we keep track of as the simulation proceeds. We let p t ,k 
represent the probability of direction reversal based on 
chemokinesis alone derived from equation (|12[l . This 
probability is then stretched towards or 1 depending 
on the local normalized chemical gradient in chemotaxis, 
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FIG. 8: Comparison of mean-field theory results with Monte- 
Carlo simulations for short times. The mean-field curves are 
evaluated after n = 600 time steps of At — 0.05 and the 
Monte-Carlo curves are determined after t — 30 time steps. 
In the chemokinesis Monte-Carlo curve $o = 4, whereas in 
the chemokinesis mean-field curve cj>o = 0.2. The factor of 
0.05 is needed to take into account the discretization of the 
numerical Volterra-type evaluation. 

through the quantity [i t ,c- Analogous to Eqs. ljT21) and 
(|26|l in the continuum theory, the kinesis and chemotaxis 
effects are represented by: 



10 5 . At these time scales, the general trends found 
in the analytic case persist for most combinations of 
(a, dp). The simulation results are summarized in a 
phase-diagram. Positive values of Sp and er both cause 
attractiveness in the system and we expect suppressed 
spreading of the distribution function. Negative values 
of dp and a on the other hand cause repulsion and en- 
hanced spreading. To quantify these effects, we compared 
the distribution function P{x, t max ) with a normalized 
Gaussian, corresponding to the parameters a = dp = 0. 
More specifically, we compare the integral of the reference 
Gaussian at the curvature inflexion points to the integral 
of the probability distribution function between the same 
limits. If the integral of the PDF is less than that of the 
corresponding Gaussian, the dynamics is classified as be- 
ing suppressed. This definition is shown in Fig.|^Jtop). 
We explore the entire phase space Sp and a and find the 
crossover from one regime to the other. The enhanced 
or suppressed diffusive character of the PDF curves may 
evolve from that of the early time regime, but it is found 
that after t = 10 3 simulation time steps the classification 
is stable. The phase diagram is shown in Fig. [3f bottom). 
For a, Sp > 0, both kinesis and chemotaxis effects sup- 
press the diffusivity of the trajectories. Only for a < 
can the spreading dynamics be enhanced. Due to our 
choice of the functional forms of the effects (Eqs. 1 121 and 
I26[l . the dynamics are much more sensitive to a that Sp, 
nonetheless, the effects of kinesis induces an asymmetry 
in Sp in the phase diagram. The curve approaches the 
limits Sp — ±0.5 asymptotically. 



1 . 
Pt.k = 77 + Sp 







(42) 



^0 t ^ t 

exp \-a{X t - Xi_i)V^] , (43) 



where V0 t = </> t (X t + 1/2) - <f> t (X t - t/2) and Sp and a 
are weight parameters subject to the constraints: 



— 2 — — 2 an< ^ — oo < a < oo. (44) 

The parameter a scales the effect size of chemotaxis, with 
no chemotaxis corresponding to a — 0; while Sp scales 
the strength of chemokinesis, with no chemokinesis cor- 
responding to Sp = 0. At each step of the Monte-Carlo 
iteration, the total reversal probability is calculated as: 



7 t =2 



1 + 



Pt,k{^t,c + 1) - 1 



Pt,k(fH,c - 1) + 1 



(45) 



In Fig. [S] we compare the results from the Monte Carlo 
simulations and the mean-field approximation of the 
previous section, for small times, for a typical case of 
chemokinesis and chemotaxis. 

We can now analyze results of the Monte-Carlo simu- 
lations for longer times than in the analytic case: t rnax = 



VI. CONCLUSION 

In this paper we presented a mean field model and 
Monte-Carlo simulations for bacterial dynamics under 
the mechanisms of chemokinesis and chemotaxis acting 
concurrently. The bacterial motion is that of a one- 
dimensional self-interacting random walker. The fibril 
trail that governs the dynamics of the system is itself a 
dynamical quantity that depends on the past motion of 
the cells. We model the two mechanisms in terms of two 
characterizing parameters, a and Sp, which represent, re- 
spectively, the degree of chemotaxis and chemokinesis. 
The mean-field results agree with Monte-Carlo simula- 
tions in the limit of short times. For long times we find a 
phase-diagram in (a - Sp) space that separates enhanced 
or suppressed diffusion regimes. In contrast to the short- 
time numerics displayed in Figs . 12171 which show PDFs 
indicative of enhanced diffusion for Sp < 0, the long time 
phase diagram exhibits suppressed diffusion provided a is 
negative and \Sp\ is not too large. Therefore, in the long- 
time limit, kinesis modeled by the logistic saturation, 
Eq. ljT21) . results in suppressed diffusion in the absence 
of chemorepellcnt effects (er < 0). We have investigated 
the effects of increased cooperativity in the biochemical 
signaling through the /? parameter. As can be seen from 
Fig|51 increased f3 values do not significantly affect the 
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does not account for contact interactions between cells, 
rather, cells interact indirectly via trails of self-secreted 
slime. This simplification enables us to study the details 
of the effects of the history-dependent fibril concentra- 
tion. Direct cell-cell signalling interactions on a 2D ag- 
gregating colony are known to lead to the formation of 
complex structures such as propagating rippling waves 
and spiraling fruiting bodies [lOt |32| . Although the in- 
vestigation of contact interactions between cells is beyond 
the scope of this paper, our approach does not preclude 
the formation of propagating waves and associated rip- 
pling patterns in systems where fibril-mediated cell-cell 
interactions are included. 

The authors are grateful for the assistance of G. 
Lakatos in improving the implementation of our simu- 
lations and numerical computations. MD and TC were 
supported by the National Science Foundation through 
grant DMS-0206733. 
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FIG. 9: (top) Gaussian and suppressed and enhanced PDFs. 
(bottom) a — Sp phase diagram at t = 10 5 for /3 = 1,4. 
The crossover between suppressed and enhanced diffusion was 
found by Monte-Carlo simulations and determining the point 
at which the integral of the reference Gaussian between the 
curvature inflexion points is equal to that of the distribution 
functions. 



phase-diagram curve. 

The model exhibits variations in shapes of the particle 
PDF. The characteristic PDF shapes depend on whether 
long or short time dynamics are considered. 

One of the possible applications of this work is for iso- 
lated A-motility (adventurous) Myxobacteria cells whose 
dynamics is driven by self-secreted slime and does not re- 
quire the presence of neighboring cells in direct contact 
|32l l33j| . Our analysis is restricted to a one-dimensional 
system. For a single cell constrained to move along a ID 
agar track, the relevant parameters (a, 5p and (3) of our 
model can be tuned to identify the correct mechanism of 
motion. Moreover, the motion of freely moving A-type 
myxobacteria cells is locally one-dimensional. Our model 
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